Three-dimensional evolution of magnetic and velocity shear 
driven instabilities in a compressible magnetized jet 

Lapo BettarinQ 

Katholieke Universiteit Leuven, Centrum voor Plasma Astrofysica, 
Celestijnenlaan 200B, B-3001, Leuven, Belgium 
Dipartimento di Astronomia e Scienza dello Spazio, 
Universitd degli Studi di Firenze, Largo E. Fermi, 2, L50125, Firenze, Italy 

Simone Landi 

Dipartimento di Astronomia e Scienza dello Spazio, 
Universitd degli Studi di Firenze, Largo E. Fermi, 2, L50125, Firenze, Italy 

Marco Velli 

Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA-91109, USA 
Dipartimento di Astronomia e Scienza dello Spazio, 
Universitd degli Studi di Firenze, Largo E. Fermi, 2, 1-50125, Firenze, Italy 

Pasquale Londrillo 
INAF Osservatorio Astronomico di Bologna, 
via C. Ranzani 1, 1-40127 Bologna, Italy 



1 



Abstract 

The problem of three-dimensional combined magnetic and velocity shear driven instabilities of 
a compressible magnetized jet modeled with a plane neutral/current double vortex sheet in the 
framework of the resistive magnetohydrodynamics is addressed. The resulting dynamics given by 
the stream+current sheet interaction is analyzed and the effects of a variable geometry of the basic 
fields are considered. Depending on the basic asymptotic magnetic field configuration, a selection 
rule of the linear instability modes can be obtained. Hence, the system follows a two-stage path 
developing either through a fully three-dimensional dynamics with a rapid evolution of kink modes 
leading to a final turbulent state, or rather through a driving two-dimensional instability pattern 
that develops on parallel planes on which a reconnection-l-coalescence process takes place. 
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I. INTRODUCTION 



Velocity shears and strong magnetic field gradients are known to play an important role 
in the dynamics of both laboratory and astrophysical plasmas, such as solar flares [H [2], 
earth's magnetotail [3111], tokamaks [5], and, in the inner heliosphere, the interaction of the 
heliospheric current-sheet with the structure determined by the slow component of the solar 
wind on the solar equatorial plane, embedded in the fast component [HI El El E] • In fact, in 
a reference frame co-moving with the slow wind, we have a bimodal flow profile where the 
velocity is zero at the HCS, and across it the interplanetary magnetic field (IMF) changes 
sign from the Southern to the Northern solar hemispheres, regions of fast wind, resulting in 
a wake flow profile whose three-dimensional evolution is not fully understood. 

The incompressible two- and three-dimensional instability dynamics of a plane current 
vortex sheet has been largely investigated in previous works flUl [TT] and only recently 
compressibility effects have started to be accounted at least for the linear regime [T2] . 
In Dahlburg et al. [10] , the simplest two-dimensional incompressible system is analyzed in the 
linear and nonlinear regime: The basic magnetic field and velocity field are both modeled by 
a hyperbolic tangent profile with a variable ratio of the velocity to the magnetic shear width, 
6, for different values of the Alfvenic Mach number, A4a- The instability behavior of such 
systems depends strongly on these two parameters. In particular, for small values of Aia 
and large values of 6, a direct energy transfer from the basic velocity field to the perturbed 
magnetic field is observed and, in general, the instability evolution can not be considered 
simply as a mixture of the Kelvin-Helmholtz (KH, hereafter) and resistive instabilities. The 
presence of a third direction is expected to modify the linear and the nonlinear regime. 
First, it allows the ignition of all the instability modes of the system and, in particular, all 
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the modes that are orthogonal to the basic flow; secondly, it allows the system to go from 
an initially laminar state to a turbulent one by means of a three-step process (a primary 
instability, e. g. the tearing mode for the neutral sheet; a two-dimensional quasi-steady 
state; a further "secondary instability", caused by three-dimensional disturbances). Such 
processes have been observed in several hydrodjTiamic and magnetohydrodynamic (MHD) 
configurations such as vortex [T3] and neutral sheets [HI [15]. Dahlburg and Einaudi [11] 
consider an incompressible problem with a subsonic shear flow whose profile is given by a 
hyperbolic tangent function with a width equal to that one of the neutral sheet. The per- 
turbations they use consist of two-dimensional disturbances of large amplitude, so that the 
system should be near the saturation condition in the nonlinear regime, with the addition 
of small amplitude three-dimensional modes. The nonlinear evolution of those systems pro- 



clockwise rotation 




Figure 1: (Color online) Schematic reconstruction of the current /neutral double vortex-sheet model 
used in the present work. The main features of the generic configuration are shown: the direction 
of the AMF, the type of rotation for the FF magnetic field (clockwise or counterclockwise with 
respect to the rotation axis, x), and the presence (absence) of the a magnetic field component 
within the fluid flow for the FF (PB). 
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duces a complex structure wherein the three-dimensional effects dominate with respect the 
two-dimensional evolution. The onset of the secondary instability has consequences also in 
the wake/jet dynamics. During the nonlinear stage of the resistive primary instability, the 
formation of magnetically confined plasma structures (plasmoids), is accompanied by the 
acceleration/deceleration of the central portion of the wake/jet [6l [7]. Once the secondary- 
instability sets in, the free energy of the velocity shear triggers the transition towards a 
turbulent state and the acceleration/deceleration stops [B]. 

Compressive effects have important consequences on the linear regime of the plane 
current- vortex sheet. Weakly evanescent modes are observed to grow and act also at large 
distances from the shear layers. Also, as shown in recent works on the three-dimensional 
evolution of the compressive tearing instability [121 E] , the onset of the secondary instability 
is observed to be critically dependent on the nature of the equilibrium, e.g. whether it is 
force-free or pressure-balanced, with or without a guide field. 

If we consider a more complex two-dimensional system given by a jet/wake interacting 
with a current/neutral double vortex sheet in a compressible and magnetically-dominated 
situation [H], the dynamics is observed to be dependent by the relative configuration of the 
basic magnetic to the velocity field, that is, as sketched in Fig. [ij by the angle a between 
the initial jet flow (or equivalently, the wake flow) direction and the asymptotic magnetic 
field (AMF, hereafter), that is the magnetic field "far enough" from the stream center. In 
particular, when a = 7t/2 in low beta configurations of supersonic, but subalfvenic shear 
flows (so, the AMF is orthogonal to the basic jet/wake flow) a typical KH instability fully 
develops and dominates the dynamics in the non linear regime. For any other angle in 
between a = (the AMF being aligned to the jet/wake) and a < 7r/2, a linear varicose- 
resistive instability [18j is dominating leading to magnetic field islands which afterwards 
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coalesce. In the present work, we extend the previous analysis to the three-dimensional case 
both in the linear and nonlinear regime. 

In the following sections, we introduce the numerical setup, the initial conditions, the 



perturbations, and the parameters for our analysis (Sect. [ITj). In Sect III, we present the re- 
sults of the simulations for the several different configurations under consideration. We will 
show how the details of the linear instability dynamics depend on the system's large-scale 
structure, in particular on the AMF's direction, regardless of the type of equilibrium consid- 



ered and how this drives into different dynamical paths in the nonlinear regime. In Sect. IV 
we summarize the results drawing our conclusions and pointing to still open questions. 



II. EQUATIONS AND NUMERICAL SETTINGS 

We solve the set of one-fluid compressive-resistive MHD equations in Cartesian geometry. 
We define a stream-wise, or Fourier, direction (y) along which we impose periodic boundary 
conditions, a cross-stream direction (x) along which the mean flow varies and we impose 
reflecting boundary conditions, and a span-wise direction (z), corresponding to an invariance 
direction for the quantities describing our initial system and along which we impose periodic 
boundary conditions. 

In a compressible MHD system where high Mach number flows are supposed to form 
and afterwards to undergo to resistively-triggered instabilities, the numerical strategy is a 
critical choice in order to obtain the best results in tracking the physical details of the ex- 
pected dynamics. On the one hand, the understanding of dissipative processes like magnetic 
reconnection requires high-order spectral-like numerical techniques such that the physical 
magnetic diffusion can be fully appreciated against the numerical dissipation. Furthermore, 
high-order methods allow to track the energy spectrum cascade towards high wave-vectors 
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with a reasonable number of grid points. On the other hand, supersonic flows can arise 
either in the initial state of the system or as a result of the reconnection process in a low 
plasma beta regime, and they determine plasma and magnetic field discontinuities that can 
be properly treated by means of shock-capturing techniques. 

Here, we use the Eulerian Conservative High-Order code (hereafter, ECHO ) to solve 
the full set of compressible and resistive MHD equations in a conservative form within the 
Upwind Constrained Transport (UCT) framework [121 ED] which properly implements the 
magnetic field divergence-free condition. In general, ECHO solves the one-fluid equations 
for a magnetized plasma in different approximations, ranging from the special and general 
relativistic MHD frameworks [2n [22] to situations wherein the presence of the magnetic field 
diffusivity is taken into account [T7] . ECHO is able to handle different high-order numerical 
techniques for the flux reconstruction and the computation of the derivatives [171 ES]. In 
particular, it allows the use of compact (or implicit) algorithms [23] which have spectral-like 
resolution properties better than the corresponding explicit numerical techniques |17j . 

The set of equations to be solved is the following: 

d,p = -V-{pw) (1) 

9t (pv) = -V- (pvv + PI-BB) (2) 
dte = - V ■ [v (e + P) - B (v • B) + 

-77 B X V X B] (3) 

dtB = - V X E (4) 

E = -vxB + r/VxB (5) 

being p the mass density, v the fluid velocity, e = p/(7 — 1) + p |vp/2 + |Bp/2 the total 
(= hydro plus magnetic) energy density according to a perfect gas equation of state, where 



p is the plasma pressure and 7 the adiabatic index; V = p + |Bp/2 the total pressure. In 
Eqs. ([3]) and (|5]), an explicit dissipative term due to the plasma resistivity, 77, is present 
though maintaining the same formal structure of the conservative framework. 

In order to obtain the dimensionless equations, we use the characteristic quantities L, p, 
Va, which correspond respectively to the velocity shear width, to the initial uniform density 
and to the the Alfven velocity at the cross-stream boundaries. Time, velocity and magnetic 
field strength are expressed in units of the related quantities i, v, B 

i=^, V = MaVa B = Va\/^T^ p. 
Va 

being M.a the Alfvenic Mach number, as already mentioned in the introduction. Moreover, 
plasma pressure and magnetic diffusivity are measured in terms of p = (3 B'^/2 and of the 
inverse of the Lundquist number S = VaL/rj. 

A. Simulation set-up 

In the present work, two different configurations are considered: a force-free (FF, here- 
after) and a pressure balance (PB, hereafter) initial equilibrium. As shown in Fig. [l| the 
FF current sheet is formed by a rotation of the cross-sheet component of the magnetic field 
across the flow, whereas in the PB case no rotation is considered and the polarity reversal is 
obtained by means of a neutral sheet. So, FF and PB cases differ essentially due to presence 
of the rotation component of the magnetic field along the flow (defined as the stream-wise 
direction) which is absent for a = (being directed along the orthogonal or span-wise direc- 
tion) and becomes more and more relevant as the angle increases. Furthermore, for cr = 
and 7r/2, the FF is intrinsically symmetric with respect to the stream- wise direction, while 
this system does not have a defined symmetry at all for < o" < 7r/2. For the PB configura- 
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tion, the system is always symmetric and in correspondence of the jet/wake we always have 
a neutral line. This initial parity property of the two configurations determines the evolution 
of the current and neutral sheet, since in the FF we have a preferential side for magnetic 
islands formation and a differential deceleration of the jet: the maximum distortion effect is 
observed for a = Svr/S 

The FF basic fields are the following 

Voy{x) = Aia sech^ (x) (6) 
Boy{x) = [sin a sech {6x) + cos a tanh (Sx)] (7) 
Boz{x) = [— cosa sech{6x) + sina tanh . (8) 

As already pointed out, current-stream interaction driven instabilities may have an impor- 
tant role in the dynamical evolution of several astrophysical plasma structures characterized 
by a low (3 regime. For instance, the slow wind acceleration region above Sun's helmet 
streamers is characterized by a typical Alfven speed less than 750 kms~^ [21], a sound 
speed of about 100 kms~^ and a typical differential velocity of fast and slow streams of 
about 300 — 400 kms~^. For our simulations, we consider the following general settings: the 
Alfvenic Mach is ~ 0.73, the sonic Mach number is equal to 3 and so we have a supersonic 
and subalfvenic current+flow system characterized by a /3 of 0.07. As already pointed out, 
the angle a defines the initial direction of the AMF relative to the basic flow, having o" = 
when the AMF is parallel to the basic velocity field and a = 7r/2 when it is orthogonal to 
the velocity field. The width of the jet/wake, a^, provides the reference length to set the 
MHD equations dimensionless, L, and hence the dynamic time is t = a^/va, while 6 = a^/aj, 
measures the ratio between the jet/wake shear width and the current sheet width (a;,). We 
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consider the perfect gas equation of state po = po Tq that initially gives 

Po = constant = 1.0 (9) 
To = = squared sound speed , (10) 

and we assume a polytropic equation with 7 equal to 5/3. 
The PB basic fields are given by the following relations 

voy{x)=MaSedi'^ (x) (11) 
Boy{x)= cos (7 tanh (5x) (12) 
Boz{x) =sma tanh (5x) , (13) 

Pressure equilibrium condition is here assured by a gradient in the temperature profile 



To = - (1 + /^oo) - - {Bl + Bl) (14) 



where Poo is the plasma beta at the cross-stream boundary for t = 0. 

The value of 6 is observed to be a critical parameter for the evolution of the above-defined 
systems and we assume it equal to 10. Such quite high value determines a strong gradient in 
the components of the basic magnetic field and this produces the conditions for an effective 
instability-triggering mechanism [Hllin])- Besides, from Eqs. ([6]) - ([s]) and Eqs. (11) - (13), 



we can observe that varying the angle a produces a rotation in the AMP. So, a strong initial 
magnetic field gradients can affect in a non-trivial way the evolution of the system according 
to the chosen initial configuration, both in the direction determined by the jet/wake and in 
the orthogonal one. 

We use two different sets of initial perturbations consisting in a proper two-dimensional 

space forcing on the cross-stream component of velocity in the jet/wake plane 

, , ^/ N . /27r/ 27rm \ 
vi^{x, y,z) = e J^[x) sm y-j- y + -j— z + <pr j , (15) 
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being (/, m) the excited modes in y and z respectively, with / G [0, 12] and m G [—6,6], (pr 
a random phase, and e the initial perturbation amplitude. In ( [l5| , we set 

= 2 tanh((5x) sech(5x) (16) 

or 

T{x) = sech{6x) (17) 

in order to force the system with a function whose parity properties should resemble that of 
the dominating instability. 

We use a three-dimensional uniform grid with resolution [n^,ny] ® = (480,120,30) 
with the following dimensionless sizes 

L, = P^[-27r,27r] (18) 
Ly = k-l [0, 2n] (19) 
L, = V,[0,2n], (20) 

where and Vz are coefficients chosen to properly adapt the size of the grid along the x 
and z directions, respectively, while kmm is the smallest perturbation wave-number along 
the ?/-direction which set the largest wavelength we can consider. The values of and 
have been chosen in order to have the cross-stream boundaries located sufficiently far from 
the current sheet in order to prevent a stabilizing effects from the walls, but at the same 
time their distance must still allows us to resolve the field gradients within the current sheet 
with a reasonable number of grid points. A reasonable value of A^^ is chosen to follow the 
dynamics both in the linear and in the nonlinear regime. Higher-resolution test simulations 
(up to double Nz) did not show sensitive differences both in the spectral behavior as well as 
in the real-space dynamics. For what regards the stream-wise direction, it must be chosen 
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Figure 2: Two-dimensional dispersion relations obtained by a linear two-dimensional code for some 
of the PB and FF cases under investigations. For the FF "7r/2" case, the system was perturbed 
by the even function (17) instead of the odd function (16). 



such to allow us to observe in the nonlinear regime both the direct cascade towards higher 
wave- vectors and a inverse cascade towards large scale by means of a coalescence process. 

We verified our choice of the main parameters {Ma, 6, (3) hj means of a two-dimensional 
linear code whose details are reported in [T7j: the maximum linear growth rates correspond to 
wave- vectors in the range [1, 3] (see Fig. [2]). We choose hence the following set of parameters 
for the numerical domain: 



6=10 



= 0.75^ = [-37r/2,37r/2] 
krmn=0.5 ^ Ly = [0,47r] (21) 
V, =1 ^ = [0, 277] . 
So, this stream- wise size implies that the wave-vector corresponding to the above mentioned 
maximum linear growth rate should corresponds to a wave-number / in a range of [2,6]. 

The simulations are performed by means of a 7^^ order compact scheme in order to have 
a low numerical dissipation influencing the expected resistivity-driven linear dynamics, and 
we consider a Lundquist number equal to 2000. As in Landi et al. [17], we verified that the 
intrinsic numerical dissipation, implicitly introduced by the adopted numerical scheme, is 
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much lower than the exphcit diffusivity assumed. In particular, by setting Uz = 1 and can- 
celing out the equilibrium magnetic field diffusion, we have verified that the two-dimensional 
linear evolution in ECHO is in very good agreement with the simulations performed by using 
the linear code. 

In the next section we present results of the three-dimensional instability both in the 



linear (§ III A) and non linear (§ IIIB) regime. To study the linear regime, as for the two- 
dimensional simulations, the equilibrium magnetic field diffusion, that affects the growth 
rates of the linear modes [I^, is cancelled out. On the contrary, the study of the non 
linear regime has been performed by including the effect of the diffusion of the equilibrium 
magnetic field. In Tab. |T| we report the simulations analyzed in the paper. 



III. RESULTS ON THE 3D INSTABILITY 
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Figure 3: Maximum growth rate as function of a for the FF (solid hne) and PB (dashed) configu- 
rations, during the hnear regime of the three-dimensional simulations. 



We consider several values of a for both the FF and PB configuration, also for continuity 
reasons with the two-dimensional analysis [U] of these configurations: 0, tt/S, ±7r/4, Stt/S, 
and 7r/2. A particular attention must be given to the FF cases with a = Svr/S in corre- 
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Table I: Summary of the performed three-dimensional simulations. In the first column, we point out 
the the two possible basic configurations: PB = "pressure balance" equilibrium; FF = "force-free" 
equilibrium. In the second column, we report considered values of the angle a. For all simulations 
the following parameters are used: 5 = 10,/5 ~ 0.07, = 3, Ada ~ 0.7,5 = 2000. 

basic configuration a 

PB, FF 

PB, FF tt/8 

PB, FF ±7r/4 

PB, FF 37r/8 

PB, FF tt/2 

spondence of which it is observed a peculiar two-dimensional dynamics characterized by a 
maximum effect in the asymmetries introduced by magnetic field basic configuration [H]. 

A. Linear regime 

Consider the results from a two-dimensional analysis: For a < 7i/2 the varicose-resistive 
modes dominate the instability dynamics [10]; a maximum growth rate of about F ~ 0.35 is 
attained for a = decreasing as a increases [9]. For a = 7r/2, the varicose-resistive modes 
are no longer unstable: the AMP is now orthogonal to the jet/wake and hence it does not 
prevent the development of a KH dynamics [HI 125] • In Fig- [2| we show the dispersion relation 
(F as function of ky) for a few cases. The a = 7i/2 case refers to the sinuous modes given 



by Eq. (17) and it is evident that the growth rate of the KH instability (F ~ 0.07) is much 



lower than those for the varicose-resistive modes with a < 7t/2. 
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Figure 4: (Color online) Magnetic energy (logarithmic scale) in the ky — Fourier space for 
three different values of a (0 in the left panel, Svr/S in the mid panel, and tt/2 in the right panel, 
respectively) for FF (upper row) and PB (lower row) at the end of the linear regime. The point- 
dashed line highlights the direction of the AMF. 

In the three-dimensional simulations, the situation changes significantly. In fact, as shown 
in Fig. [3| we obtain the same growth rates of about F ^ 0.35, regardless of the value of 
a. Moreover, all the cases are characterized by the development of varicose-resistive modes, 
as we can infer from the fact that the growth rates match the values observed in the two- 
dimensional linear simulations with odd perturbations. In general, the duration of the linear 
phase is slightly longer in the PB cases than in FF's ones since the growth rates of the fastest 
modes are a bit smaller. This is consistent with the 2D linear simulations. This different 
behavior is clearly related by the fact that in three dimensions all instability modes are 
allowed to developed. In Fig. |4| we report the magnetic energy spectrum in the Fourier 
a given time during the linear regime for the FF (upper row) and the PB 
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Figure 5: From left to right, the most relevant Fourier modes of the magnetic energy in the nonlinear 
phase as a function of time for the FF configurations with a = (left), a = tt/A (middle), and 
0- = 7r/2 (right). 

(lower row) configurations: the cases with cr = 0, Svr/S and 7r/2 are shown in the left, middle, 
and right column, respectively. As the instability dynamics sets in, both the FF and the 
PB configurations where the AMF is aligned with the basic fluid jet (a = 0) are driven 
mainly by the modes characterized by = 0, although a small amount of energy is present 
in modes with kz = ±1, 2, at least for the PB configuration. Consistently with the linear 
code results. The observed fastest-growing linear mode lies between ky = 2, 2.5. If now we 
consider a value of cr between and tc/2, such as a = Sir/ 8 shown in the middle panel of 
Fig. |4| we observe that the most unstable modes are no more aligned with the jet/wake 
flow. The instability grows along a preferential direction that is selected by the direction of 
the AMF, underlined in the figure by the dash-dotted line. The most unstable modes for 
this case are those characterized hj ky = 1 and k^ = 2,3, corresponding to \k\ in the range 
[V^, "\/T0] which is consistent with the values predicted by the linear simulations for a = 0. 
Same results, not shown here, are obtained for a = 7r/8, and a = n/A. Moreover, the same 
results is obtained in the case of a = — 7r/4 (FF) configurations: this implies that the way 
the cross component of the magnetic field rotates through the current sheet is not relevant 
for the dynamics of the system. 
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As shown previously, in the two-dimensional case, for a = 7r/2, the system should be 
unstable under KH like instability. However, if we look at the right panel of Fig. |4| we note 
that both FF and PE configurations are dominated by modes aligned with the z-direction, 
the direction of the AMF field, although a significant dose of energy is still present in 
modes with kz = for the PB case. The most unstable modes lie in the wave-vector range 
kz ~ ±2 — 3: the system is dominated by resistive-varicose instabilities. The presence of 
excited modes along the y-direction is due to the presence of the KH instability, stronger for 
the PB case because of the absence of the stabilizing effect of the perpendicular magnetic 
field inside the current-sheet. 

In conclusion, regardless of the a values, the most unstable modes driving the system 
throughout the linear regime are resistive-varicose and they are selected according to the 
AMF direction regardless of the detailed structure of the magnetic field within the current- 
sheet, which is the type of equilibrium under consideration. In particular, these modes are 
characterized by the condition /c x B(|5x| ^ 1) = 0. Such preferential direction determined 
by the asymptotic magnetic field in the linear evolution of the instability was observed also 
in the case of the incompressible and compressible three-dimensional tearing mode evolution 
by Onofri et al. [I6] and Landi et al. [I7j respectively. 

B. Nonlinear regime 

Although the nonlinear regime of the current (neutral) double vortex sheet is charac- 
terized by a complex behavior which depends on the initial evolution of the most unstable 
modes, it is possible to identify some key-features. 
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Figure 6: (Color online) Magnetic energy in the {ky — k^) Fourier space for the FF configurations 
corresponding to different values of a in the well-developed nonlinear regime. The dashed line 
highlights the o" — 7r/2 direction. 



Figure 7: (Color online) Three-dimensional plasma pressure iso-contour at t = 200 for the FF 
configuration with a = 0. The shaded surface encodes the region of space where the pressure is 
more than 65% of its maximum. 
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Figure 8: (Color online) Three-dimensional plasma pressure iso-contour at three different instants 
during the nonlinear evolution of the FF configuration with a = 7r/4. These images refer to t = 100 
(left panel), t = 150 (middle panel) and t = 200 (right panel). The shaded surface encodes the 
region of space where the pressure is more than 50% of its maximum. 

1. FF cases 

We start from Fig. |5] to consider all the cases determined by a different value of a. 
Here, we show the temporal evolution of some magnetic energy modes for three different FF 
configurations: a = (left panel), a = ti/A (middle panel) and a = tx /2 (right panel). 

a = 0: From about t ~ 50, it is clearly observed a coalescence process driving the system 
to the maximum length scale allowed by our numerical domain. In fact, the inverse 
cascade ends as soon as the largest wave-number mode dominates at about t ~ 80. 
In the meanwhile, the orthogonal modes (those with fcj^ = and kz 7^ 0) keeps on 
growing and, from t ~ 150 on, they start to determine the structure of the system. 
At about t ~ 200, most of the energy is symmetrically distributed in the "secondary" 
(orthogonal) modes as well as in the first primary modes {ky 7^ and kz = 0), as 
shown in the top-left panel of Fig. |6] where the magnetic energy spectrum in the plane 
{ky, kz) in the well-developed nonlinear regime is shown for different a-valued FF initial 
configurations. This dynamics is analogous to those shown in Landi et al. [IT] analyzing 
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Figure 9: (Color online) Three-dimensional plasma pressure iso-contour at the end of the simulation 
(t = 200) for two different FF configurations: a = tt/8 (left panel) and a = Svr/S (right panel). The 
shaded surface encodes the region of space where the pressure is more than 55% of its maximum. 



Figure 10: (Color online) Three-dimensional plasma pressure iso-contours of the FF configuration 
with cr = tt/2. Image refers to t = 125 (left panel) and t = 200 (right panel). The shaded surface 
encodes the region of space where the pressure is more than 50% of its maximum. 

the three-dimensional compressive tearing evolution. Also in that study, the nonlinear 
magnetic energy spectrum is characterized by two preferred directions resulting in a 
system modulated in two different ways in the physical space. In Fig. [7] it is shown 
the three-dimensional plasma pressure iso-contour for the "0" case at t = 200 and 
it is possible to observe a coalesced plasma structure along the stream-wise direction 
modulated along the span- wise axis. 

< a < 7t/2: The overall dynamics is not very different from the previous case, provided 
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Figure 11: From left to right, the most relevant Fourier modes of the magnetic energy in the 
nonlinear phase as a function of time for the PB configurations with a = (left), a = 7r/4 
(middle), and a = ix jl (right). 

that the nonlinear driving modes are properly considered. For instance, in the "7r/4" 
case shown in Fig. |5] (middle panel), we observe a coalescence process with about the 
same time-scale of the "0" case. Afterwards, the modes with a given direction to the 
AMF starts to increase, in particular the mode (1, —1), and in the far nonlinear regime 
they own most of the system energy as shown in the top-middle panel of Fig. [6] In the 
physical space, the rise up of the secondary instability determines that the globally- 
coalesced structure observed in the left panel of Fig. [8] at about t ~ 100, at later 
times is modulated along the direction located at 7r/4 with respect the x — y plane as 
shown in the middle and right panels of Fig. [8j In Fig. |6] we show also the magnetic 
energy spectra at the end of the nonlinear regime also for the cases a = 7r/8 (bottom- 
left panel) and Svr/S (bottom- right panel). The dashed lines represent the cr — 7r/2 
direction and they well fit the observed spectrum anisotropy direction. In all cases. It 
is worth pointing out that this privileged direction at about t ~ 200 is orthogonal to 
the direction of the most energetic modes in the linear regime (see Fig. [2] as reference 
for the 37r/8 case). The three-dimensional plasma pressure iso-contours at the end 
of the simulations for the a = vr/S (left panel) and for a = Svr/S case (right panel) 
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shown in Fig. |9] reveal the effects in the physical space of the behavior observed in the 
magnetic energy spectra. 

a = As for the previous cases, the first part of the non linear evolution is characterized 
by the growth of the primary instability along the AMF, that is now the span-wise 
direction z. At about t = 125, the largest coalesced structure with respect to our 



numerical box in the z direction is developed. As shown in the left panel of Fig. [TO 
in the physical space enhancements of the plasma pressure underlines the presence 
of magnetic islands, as it is typical of resistive driven instabilities. Yet, the modes 
orthogonal to the primary ones are still growing and they reach the same energy 
level at about t = 150. So far, system's behavior is mostly analogous to the previ- 
ously described cases, characterized all by an energy density almost equally distributed 
both in primary and secondary modes. Yet, the final part of the nonlinear evolution 
presents a completely different dynamics: the modes along the stream-wise direction 
are dominating the others of about one order of magnitude in the energy spectrum (see 
Fig. [5] and |6]). An analysis in the physical space reveals that the prevailing modes have 



a different nature with respect the previous cases. In fact, in the right panel of Fig. [TO 
we observe that at t = 200 the plasma pressure shows the typical sinuous profile of 
KH instabilities. So, although it is completely overwhelmed by the resistive- varicose 
mode in the linear stage, nevertheless the KH dynamics appears to dominate in the 
well-developed nonlinear regime and as it seems to be mainly driven by the presence 
of the velocity shear, we expect it to influence strongly the wake (jet) acceleration 



(deceleration) (see in section IIIC). 
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So, as already mentioned, we observe a qualitatively similar behavior for different values 
of 0" < 7r/2. A resistive linear regime is driven by the modes selected according to the 
direction of the equilibrium asymptotic magnetic field initially defined in the (y, z) plane. 
The nonlinear regime is initially characterized by the formation of magnetic islands then 
undergoing to a coalescence process. Afterwards, the onset of secondary instabilities kink 
the coalesced structures: in the later stages of the simulations, we observe plasma pressure 
enhancements strongly modulated by modes orthogonal to the primary modes, as similarly 
observed in three-dimensional tearing simulations [17] . In the Fourier space, the presence of 
the secondary modes is evident in the magnetic density energy contained in modes orthogonal 
to the primary ones. For o" = 7r/2, the system is also unstable to KH instabilities and this 
changes significantly its nonlinear dynamical path. In this case, the later stage of system's 
evolution is dominated by sinuous modes. 



2. PB cases 



Now, we consider the nonlinear evolution of the PB equilibrium configurations for different 



values of a. In Fig. 11 , it is shown the time evolution of the most important magnetic energy 



modes for the cases a = 0, cr = 7r/4, and a = n/2. 

(7 = 0: We observe that at the beginning of the nonlinear regime the system is still driven 
by one of the linearly-dominating modes, {ky,kz) = (1,0) leading it to be structured 
almost on the largest length-scale. In fact, after a few dozens of time steps we observe 
a inverse cascade towards the maximum length scale allowed by our numerical box, 
(0.5,0), followed by a saturation plateau of the nonlinear regime lasting up to about 
t ~ 130. In the meantime, the orthogonal modes (0, ±1) and (0, ±2) grow as well. So, 
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the two-dimensional structured system starts to develop in the z-direction according 
to these modes which, together with their harmonics, drive the system till the end of 
the simulation. Primary modes directed along the ?/-direction, except for the funda- 
mental one (0.5,0), differently from the FF configuration, are strongly damped: This 



is evident, for example, by comparing the mode (1.0, 0) in Fig. 11, left panel, with the 
analogous one for the FF configuration in Fig.|5} So, from an early coalescing structure, 
the system presents at the end of the simulation an essentially two-dimensional con- 



figuration in the ( observed in the plasma pressure structure shown in Fig. LL2 



This behavior is confirmed also by the magnetic energy spectrum at the end of the 



simulation, shown in the left panel of Fig. 13: most of the magnetic energy is now 
contained in modes directed along the ^-direction and we do not observe the energy 
redistribution between primary and secondary modes as in the analogous FF config- 
uration. This result is consistent with those obtained in the pure three-dimensional 
tearing instability simulations for PB configuration |17j . 

< a < 7r/2: Let us focus on one emblematic case, such as for example a = tt/A. After the 
first phase, where the system evolution is determined by the primary three-dimensional 



we 



modes characterizing the linear regime, as shown in the middle panel of Fig. 11 
have an increase of the contribution of the mode {ky, k^) = (0.5, 0) which drives the 
dynamics till the end of the simulation. At the same time, modes orthogonal to the 
primary ones increase in time such that, at the end of the simulation, they contain a 
large portion of the magnetic fluctuating energy, as it is shown in the middle panel 



of Fig. 13 Again, the dashed line locates the a — n /2 direction and it well fits the 
observed spectrum anisotropy direction pointing out that the privileged direction at 
about t ~ 200 is orthogonal to the direction of the most energetic modes in the linear 
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Figure 12: (Color online) Three-dimensional plasma pressure iso-contour during the nonlinear 
evolution of the PB configuration with o" = (left), o" = 7r/4 (middle) and u = 7r/2 (right). The 
shaded surface encodes the region of space where the pressure is more than 60% of its maximum. 
For (7 = and a = 7r/4 the pressure is evaluated at t = 200, while for the a = tt/2 the reference 
time is t = 150. 




Figure 13: (Color online) Magnetic energy in the {k^ — ky) Fourier space for the PB configurations 
with different values of a at i ~ 200. As for the FF case (Fig. [6|, the dashed line highlights the 
angle a — 7r/2. 

regime. A similar behavior is observed in the magnetic energy spectrum for the "vr/S" 



and "37r/8" cases (bottom-left and bottom-right panels of Fig. 13, respectively). The 
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physical behavior of the system is shown in Fig. 12 (middle panel), where it is evident 



the instability modulation of the whole structure along the preferred direction. 

a = tt/2: The overall behavior is similar to the corresponding FF configuration, although it 
appears faster in its evolution. At t = 50, it is attained the saturation in the magnetic 
energy power in the smallest span-wise (z) wave-vector and it is observed a rapid 
growth of the energy of stream-wise (y) modes. At t = 100, the system is strongly 
structured in the {x,y) plane with a typical varicose-mode configuration. However, 
the presence of KH instability modes leads the system to be structured in a typical 
sinuous configuration already at t = 150 as it is possible to observe in the right panel 



of Fig. 12 At the end of the simulation, most of the magnetic energy is contained in 



stream-wise (y) wave- vectors (Fig. 13, top right panel). 



C. Jet /Wake evolution 

In general, a net acceleration effect along the stream-wise direction of the jet/wake em- 
bedding the current sheet is observed, while the mean (positive and negative) velocity con- 
tributions along the other directions are negligible for almost all the nonlinear regime. We 
can measure the net acceleration effect by averaging out the span-wise and stream-wise de- 
pendence of the stream-wise component of the velocity field, i. e. by looking at the function 

{vy)yz = -j—j- / Vy{x,y,z)dydz. (22) 

J^yJ^z Jo Jo 



In Fig. 14, we plot this quantity for different values of the angle a and for both FF and 
PB equilibria (first and second panel). For comparison, we report the same profiles for 
two-dimensional simulations with the same initial condition in the well-developed nonlinear 
regime (third and fourth panel). It is manifest the strong differences from the two- and the 
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Figure 14: (Color online) Wake profile of the stream- wise component of the velocity, as computed 
using Eq.|22[ for all the values of a. The black solid line in all plots is the wake profile a t = shown 
as reference. The first and second panels refer respectively to the FF and PB three-dimensional 



cases (see section IIIC), whereas, for comparison, in the third and fourth panels it is reported the 



FF and PB two-dimensional cases respectively. 

three-dimensional case. Even though we are exploring a slightly different parameter space 
from Bettarini et al. [9J, the two-dimensional simulations we performed confirms the overall 
dynamics there discussed. In three-dimensions, we observe a weaker acceleration effect and 
a less pronounced enlargement of the wake, both in the FF and PB cases. Because of 
the insurgence of the secondary instability, dips in the velocity profiles observed in two- 
dimensional simulations, caused by the presence of coherent vortices structures inside the 
current-sheet region, are no longer present. 



IV. DISCUSSIONS AND CONCLUSIONS 

We consider the fully three-dimensional linear and nonlinear evolution of a plane neu- 
tral/current double vortex sheet consisting in a neutral/current-sheet embedded in a sheared 
supersonic flow. Several parameters can influence the evolution of the initially perturbed 
system: The initial equilibrium configuration, the magnetic and velocity fields relative ge- 
ometry, the plasma /5, the sonic and Alfvenic Mach numbers. Here we have considered the 
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typical values describing jet/wake flows in the solar active regions at a few solar radii: in 
particular, the conflguration can be applied to the wake model of the slow component of the 
solar wind ^UU- 

Two different equilibrium conflgurations have been compared: a force-free (FF) magnetic 
fleld whose polarity reversal determines the formation of a current sheet inside the plasma 
flow and a pressure balance (PB) magnetic fleld where the polarity reversal is obtained 
by means of a neutral plane inside the flow. Differently from the 2D simulations [9j, for all 
orientations of the asymptotic magnetic fleld, both the FF and PB equilibrium conflgurations 
are dominated by varicose-resistive modes. The most unstable ones are characterized by 
wave- vectors whose direction is parallel to the asymptotic magnetic fleld. The presence of a 
rotational component inside the current sheet, as for the FF case, does not affect signiflcantly 
the linear behavior of the system. Such selection rules has been already observed in previous 
analysis of three-dimensional tearing mode evolution [T7] . 

The nonlinear regime is always characterized by a two-stage evolution: for a < 7r/2, 
the varicose-resistive modes (primary modes) start to coalesce and then a secondary insta- 
bility develops. In the physical space, it results in the formation of magnetic islands and 
pressure enhancement regions which afterwards are deformed by kink-like instabilities. The 
appearance of the secondary instability is accompanied by the growth of fluctuations whose 
wave- vector are orotogonal to the primary ones. The resulting magnetic energy spectrum, as 
for the simple current-sheet instability [ISIIIZI, results thus strongly anisotropic. As pointed 
out by Landi et al. fT7], the presence of the magnetic fleld inside the current sheet in the 
FF conflguration can reduce the effects of the secondary instability: as a consequence, in 
the fully developed non linear regime a large dose of magnetic energy is in both the primary 
and secondary (orthogonal) modes. Because of the absence of a guide fleld, which has a 
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stabilizing effect on the current-sheet evolution [T6t [TTf [26] , for the PB configuration, the 
transition to the secondary instability appears earlier and the magnetic energy spectrum is 
dominated by these modes. As in Bettarini et al. (H] the production of magnetic islands is 
accompanied by vortices formation which are considered the responsible of the wake (jet) 
acceleration (deceleration) [6l [7]. Similarly to what observed by Einaudi et al. [6], the on- 
set of the secondary instability determines the destabilization of the islands, and hence it 
reduces the acceleration/deceleration effect. 

For a = 71/2, for sonic and Alfvenic Mach numbers typical of the heliospheric environ- 
ment, the system results to be unstable both by tearing and Kelvin-Helmholtz instabilities. 
With the relatively low Lundquist number here adopted, it results that, differently from 
the two-dimensional case [Hj, in three dimensions, the linear regime is dominated by resis- 
tive instabilities. However, during the non linear regime, because of the dissipation of the 
underlying equilibrium magnetic field, the width of the current-sheet increases and, as con- 
sequence, there is a decrease in time of the tearing mode growth rate as it goes as 5^^"^ |27] . 
As a consequence, if the diffusion acts for a sufficiently long time, the flow driven instability 
is able to dominate the system during the non linear regime. Moreover, since the growth rate 
of the resistive mode scales as iS~^/^, it is expected that, for a magnetic configuration where 
the asymptotic magnetic field is orthogonal to the jet/wake, Kelvin-Helmholtz instabilities 
will be the prominent dynamics driving the overall system evolution. In the present work 
the current (neutral) double vortex sheet evolution is investigated in the limit of zero kine- 
matic viscosity, in analogy with [TTj. In fact, although the influence of viscous effects can be 
retained with respect to the resistivity-driven dynamics, however it has been suggested that 
the overall evolution depends essentially on the Hartmann number which gives a measure 
of the relative importance of drag forces resulting from magnetic diffusivity and kinematic 
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Figure 15: (Color online) Top panel: Volume average of the ohmic and viscous dissipation (solid 
and dashed line, respectively) as a function of the simulation time for the FF simulation with 
a = 7r/4. Bottom panel: Ratio of these quantities as a function of the simulation time. 



viscous forces [28]. Fig. 15 shows an a posteriori computation of the volumetric average of 
the ohmic and the viscous dissipation (sohd and dashed hne respectively, in the top panel) 
and their ratio as a function of the simulation time. Here, a FF simulation with cr = 7r/4 is 
considered, but this case is a paradigma for all simulations. The above mentioned quantities 
are defined respectively as 



y\^oh^{t)= / dV7]\3Y (23) 
Jv 

^^v^s{t)= [ dvu\n\\ (24) 



Jv 

where we assume u = rj, Q is the vorticity vector, V x v, and V is the volume of our three- 
dimensional numerical box. In general, the overall ohmic dissipation is far bigger than the 
viscous one: as the motion gets turbulent and vortices form and grow, the contribution of the 
viscosity becomes more and more important up to be one fifth of the other. Yet, the overall 
dynamics throughout the linear and nonlinear regime is not affected by that, in particular 
with regards to the onset of the coalescence process or the triggering of ideal instabilities 
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as the secondary ones [TT]. Nevertheless, a further investigation on viscosity effects in 
longer simulations is mandatory. For several heliospheric and astrophysical environment 
applications it will be suitable to consider the effect of the radial expansion and/or the 
presence of an underlying stratified medium: they can influence not only the onset of the 
primary instability, but also the length- and time-scales associated with the insurgence of 
the secondary instability. 
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